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A general lattice Boltzmann (LB) model is proposed for solving nonlinear partial differential equa- 
tions with the form dt<fi + XX=i a kd x Ilk{(j>) = 0, where at are constant coefficients, and Hk((fr) are 
the known differential functions of <f>, 1 < k < m < 6. The model can be applied to the common 
nonlinear evolutionary equations, such as (m)KdV equation, KdV-Burgers equation, K(m, n) equa- 
tion, Kuramoto-Sivashinsky equation, and Kawahara equation, etc. Unlike the existing LB models, 
the correct constraints on moments of equilibrium distribution function in the proposed model are 
given by choosing suitable auxiliary-moments, and how to exactly recover the macroscopic equa- 
tions through Chapman-Enskog expansion is discussed in this paper. Detailed simulations of these 
equations are performed, and it is found that the numerical results agree well with the analytical 
solutions and the numerical solutions reported in previous studies. 

PACS numbers: 02.70.-c, 02.60.Cb, 05.45.Yv 



I. INTRODUCTION 

The lattice Boltzmann method (LBM) is a promising technique for simulating fluid flows and modeling complex 
physics in fluids [J, HI, 9 • Compared with the conventional computational fluid dynamics approaches, the LBM is easy 
for programming, intrinsically parallel, and it is also easy to include complicated boundary conditions such as those in 
porous media. Up to now, the most widely used LBM is the so-called lattice Bhatnagar-Gross-Krook (LBGK) model. 
However, the LBGK model may suffer from numerical instability when it is used to simulate the fluid with small 
viscosity. A lot of work has been done to improve the stability of LB model, among which the multi-relaxation-time 
LBM $MM and entropic LBM 0, [1, i, [M 

1 1 XI ] have attracted much attention in recent years. It should be noted that 
the LBM also shows potentials to simulate the nonlinear systems, such as the reaction-diffusion equation fl2l. [T3L fl3| . 
convection-diffusion equation (CDE) [H, [H, [13, Ell j Burgers equation [b|, KdV-like equation [2(|, Poisson equation 
[2ll ] . etc. Recently, the LB models have been extended to solve CDEs on rectangular or irregular lattices [I2,[23[ and 
anisotropic dispersion equations [UHf! Hsj], among which the model proposed by Ginzburg [26[ is generic. 

Except for solving real- valued nonlinear systems, the LB and LB-like models have been successful in solving complex- 
valued nonlinear systems. Since the middle of 1990s, several types of quantum lattice gases and quantum LBM have 
been proposed based on quantum-computing ideas to model some real and complex mathematical-physical eq uations, 
such as the Schrodinger equation, Gross-Pitaevskii equation, Burgers equation, KdV equation [27], [H, [H, Ho, HH, [H, 
[H, [H, [H, HH, etc. We refer the readers to a recent paper [3(| for a detailed review. On the other hand, recently 
the classical LB model has been used to model complex-valued equations. In Ref. [37| the LBM was applied to 
one-dimensional (ID) nonlinear Schrodinger equation (NLSE) following the idea of quantum lattice-gas model (30l.[3l| 
to treat the reaction term. In Ref. (3a |. motivated by the work in Ref. [37], the LBM for n-dimensional (nD) CDE 
with a source term was directly applied to some nonlinear complex equations, including the NLSE, coupled NLSEs, 
Klein-Gordon equation and coupled Klein-Gordon-Schrodinger equations, by adopting a complex-valued distribution 
function and relaxation time. In Ref. [39(, a general LB model for a class of nD nonlinear CDEs was presented by 
properly selecting equilibrium distribution function. The model in Ref. [39] ] can be applied to both real and complex- 
valued nonlinear evolutionary equations. Following the idea in Ref. (39j . a LB model for ID nonlinear Dirac equation 
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was given in Ref. (4o| . which is of second-order accuracy in both space and time, and the order of accuracy is near 
3.0 at lower grid resolution. The studies in Refs. [13, [H, H^, H(| show that the LBM may be an effective numerical 
solver for real and complex- valued nonlinear systems. 

Most of the existing LB models are used for solving partial differential equations (PDEs) with order lower than or 
equal to three, while many efficient conventional numerical methods for solving higher-order PDEs have been proposed, 
such as the pseudo-spectral method local discontinuous Galerkin (LDG) method 0, S|, finite different scheme 
[44| . finite volume method (4f|, radial basis function method 46], among them LDG method has attracted much 
attention due to its good nature, such as flexibility and high parallel efficiency (see a recent review article [43| and 
references therein for details) . Although some higher-order LB schemes have been proposed, they are mainly limited 
to solve lower-order PDEs [47],[48|, [4f| and ID special problems [50]. Furthermore, the efficient numerical analysis of 
these schemes are still needed. Therefore, it is important to research and develop the LB model for solving higher-order 
PDEs. 

In this paper, by extending the idea in Ref. [39[ a general LBGK model is proposed for solving a class of nonlinear 
partial differential equations (NPDEs) with order up to six. In order to exactly recover the macroscopic NPDE, the 
correct constraints on moments of equilibrium distribution function in the proposed model are given by introducing 
suitable auxiliary-moments. The proposed model can be used to solve many common nonlinear evolutionary equations, 
such as (m)KdV equation, KdV-Burgers equation, K(m,n) equation, Kuramoto-Sivashinsky equation, and Kawahara 
equation, etc. Numerical results show that the LBGK model can also be used to simulate higher-order NPDEs. 

The rest of the paper is organized as follows. In Sec. II, the LBGK model for NPDE is presented, and how to 
exactly recover the macroscopic equation from the model is discussed. In Sec. Ill, the equilibrium distribution and 
auxiliary-moment functions for NPDEs with different orders are given. Numerical tests of the LBGK model are made 
in Sec. IV, and finally a brief summary is given in Sec. V. 

II. MULTI-SCALE LATTICE BOLTZMANN EQUATIONS 
A. LBGK Model 

The ID NPDE considered in this paper can be written as 

m 

d t cj> + Y J O t kd k Jl k {cj>)=Q, (1) 

k=l 

where is a scalar function of position x and time t, a\ = 1, a/j are constant coefficients, and Hk{4>) are known 
differential functions of <^>, I < k < m < 6. 

Our LBGK model is based on the D1Q& lattice [2] with b velocity directions in ID space. The evolution equation 
of the distribution function in the model reads 

fj{x + Cj At, t + At) = fi(x, t) - -[fjix, t) - /;>, t)], (2) 

r J 

where {cj, j = 0, ...,b— 1} C {0, c, — c, 2c, — 2c, 3c, — 3c} is the set of discrete velocity directions, c = Ax /At, Ax 
and At are lattice spacing and time step, respectively, r is the dimensionless relaxation time, and /J q (x,t) is the 
equilibrium distribution function (EDF). 

To solve Eq. §Q using Eq. we must give appropriate EDF fj q (x,t). By applying the idea in Ref. [39( to the 
higher-order NPDE (TTJ), the following constrains on /J q are given 

Y / ^fP = n k0 +(3 k Il k ,k = 2 1 ... 1 m (3) 

j 

where I3 k are parameters, and n^o are auxiliary-moment (AM) functions for correctly recovering Eq. {l}, which are 
determined later, k = 2, . . . , m. 
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B. Multi-Scale Lattice Boltzmann Equations 



To derive the macroscopic equation ([I]), the Chapman-Enskog (C-E) expansion in time and space is applied: 

n = E ek ti k) * dt = E ^ - 5 * = e ^ > ( 4 ) 

fe=0 fe=l 

where e is a small expansion parameter. Using the first formula in Eq. ([3]) and Eq. [[I"]), we have 

E/f )= M>i. (5) 

3 

By applying Taylor expansion to Eq. ([2]), we get 

D & + f ^/3 + • • • + + ■ = ~ifi - (6) 

where Dj = dt + Cjd x . Denote D\j — dt x + Cjd Xl . Similar to Ref. [13], substituting Eq. Q into Eq. {B} and treating 
the terms in order of e k separately gives 







if = f-\ 


(7) 






1 f (i) 
rAt Jj ' 


(8) 


d t2 ff 


+ r 2 A«^/f = 


1 f (2) 

rAt^' ' 


(9) 


d t jf ] +2r 2 Atd t2 D lj j fU 




1 f (3) 

rAt Jj ' 


(10) 


d t jf ] + 2r 2 Atd t3 D l3 ff + 3r 3 At 2 d t2 DyP + n&Kf^ + 


-r 4 At 3 Dy^ = 


1 f (4) 


(11) 



d t jf + 2r 2 Atd u D x ,!Y> + frnjAi^By/^ + ^At'd^ff 
+2r 2 Atd t2 d t jf ] + Ar^dt^ff + r 5 At 4 I^./j 0) = -^/f, (12) 
d t6 ff + 2r 2 Ata t5 i? lj /j 0) + 6r 3 At 2 9 t2 a t3 i? lj /j 0) + ^d^dul^ff + 2r 2 Aid t2 d t4 /j 0) + 6r 4 A* 3 ^ 2 D i ! J ./j 0) 

+r 3 Ai 2 a t 3 2 /j 0) + r 2 A^ 3 /j 0) + 4r 4 At 3 9 t3 ^/f + 5r 5 Ai 4 3 t2 ^/j 0) + r 6 At 5 Dyf = (13) 



where 



T 2 = 


- r4 


" 2' 




T 3 = 


T 2 


1 




T 4 = 


- T 3 


3 2 7 

+ — T T 4 

2 12 


1 

" 24' 


?5 = 


T 4 


- 2r 3 + ^r 2 - 
4 


1 

4 T + 


T6 = 


- T 5 


, 5 4 13 3 
H T T 


3 2 




2 6 


+ r 



1 

120' 

-^ + J_. (14) 
360 720 v ' 



C. Recovery of the Third-Order NPDE 

Since general LB models for the second-order NPDE have been developed (2fj| . [3^ | , we only discuss how to recover 
NPDEs with orders higher than two in this paper. Summing Eqs. ([5]) and © over j, and using Eqs. ([3]) and (JSJ), we 
obtain 

E Dijff =d tl + d^U^) = 0, (15) 

j 

d^ + ^At^Dyf -0. (16) 
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Using Eqs. ((3]) and (fTB")) , and taking n 2 o as in Ref. [39| such that <9 tl IIi + d Xl Tl2o = 0, that is n 2 o = / djjlidjjlidcf), 
we have 

J2 D ljf? ] = d l ( t ) + 2a tl 5 Xl ni + dl x (n 20 + /3 2 n 2 ) - 9 Kl (a tl n x + s xi n 20 ) + ftfi^ifc = ft^n 2 , (17) 
j 

then it follows from Eqs. (HI and (JT7J) that 

St^ + aa^na = 0, (18) 

with a 2 = Air 2 /3 2 . 

Summing Eq. (fTU|) over j, and using Eqs. (5]), (US]) and (fT7|) . we obtain 

9 t3 + r 3 At 2 ^^/jo )=O; (1Q) 

i 

]T D L /i 0) = 9 l <t> + 3d l d *i n i + 3 ^ C ( n 2o + &n 2 ) + ^ (n 30 + /3 3 n 3 ) 

= 2^^^ + 3d tl d 2 xi (n 20 + /? 2 n 2 ) + d 3 xi (n 30 + /3 3 n 3 ) 

- ^ 1 (9 tl (n 20 + 3/3 2 n 2 ) + 9 :ci n 30 ) + /3 3 93 i n 3 . (20) 

If we take IT 3 o such that c* tl (n 20 + 3/3 2 II 2 ) + d Xl U 30 = 0, that is, 

n 30 = J d^(u 20 + 3p 2 Ti2)d^u l d4> = J [(^n x ) 3 + 3/3 2 5 n 2 a n x ] (21) 

then from Eqs. (fH?|) and (j2"0|) , we have 

c^ + a^II^O, ( 22 ) 

with a 3 = Ai 2 r 3 /3 3 . 

Combining Eqs. (jT5j> , ([15]) and (|2"2"|) , the third-order NPDE is exactly recovered to order 0(e 3 ) 

9 t + ^n 1 (0)+a 2 ^n 2 (0)+a 3 93n 3 (0) =0, (23) 

with a 2 = Air 2 /3 2 , a 3 = Ai 2 r 3 /3 3 . 

Remark 1. Note that there are no additional assumptions on the present model for the third-order NPDEs with 
the form of Eq. (Q]) , and if we select II 2 o and II 3 o above the third-order NPDEs are exactly recovered. The present 
model with a D1Q4 or D1Q5 lattice can be used to simulate the third-order NPDE ((T|) which contains some KdV-type 
equations. 

D. Recovery of the Fourth-Order NPDE 
Summing Eq. (TT]) over j, and using Eqs. ©, ©, (US]), (JTTJ) and we obtain 

9 t4 + S^Ai 2 ^^ (/3 2 n 2 ) + r 2 Aid 2 2 + r 4 At 3 Dyf ] = 0. (24) 

3 

Y, D hfP = %d> + wldvji! + 6dld 2 xi (n 20 + /3 2 n 2 ) + Ad tl d 3 xi (n 30 + /3 3 n 3 ) + d 4 xi (n 40 + /3 4 n 4 ) 

3 

= zdldvjix + 6dld 2 xi (n 20 + /3 2 n 2 ) + 4d tl d 3 xi (n 30 + /3 3 n 3 ) + s£ (n 40 + /3 4 n 4 ) 
= (3n 20 + 6/? 2 n 2 ) + 4^^ (n 30 + /? 3 n 3 ) + ^ ( n 40 + /? 4 n 4 ) 
- 9^,^20 + d tl d 3 xi (2n 30 + 4/3 3 n 3 ) + fl£ (n 40 + /? 4 n 4 ) 

= fig, (9 tl (2n 30 - n 30 + 4/? 3 n 3 ) + a Kl n 40 ) + /^n,, (25) 
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where II 3 o satisfies that d ti n 20 + d Xl Il3o = 0, that is, IT30 = J d^Uwd^Uidcf). 
If we take IT40 such that 9 tl (2n 30 — II30 + 4/3 3 II 3 ) 4- fl^ILio = 0, that is, 

n 40 - fd^Uao - 1T30 + 4 / 9 3 n 3 )^n 1 # = J [(c^il) 4 + 6p 2 d 4> n 2 {d 4> n 1 ) 2 + i^d^d^] d<j>, (26) 

then from Eqs. (|24|) and (j25|) we have 

d ti 4> + 3T 3 fcAt 2 d t2 dljl 2 + T 2 Atd? 2 + and^JU = 0, (27) 

with «4 = At 3 T4/34. 

Note that there are two additional terms 3r 3 [3 2 At 2 dt 2 d 2 Jl 2 and T 2 Atd 2 2 (f> in Eq. (|27[) should be removed when 
correctly recovering Eq. {T]). From Eq. (fT8|) . we have 

3T 3 [3 2 At 2 d t2 d 2 Xl Tl 2 + r 2 Atd 2 2 4> = (3/3 2 r 3 At 2 - a 2 T 2 At)d t2 d 2 Xl U 2 = (3r 3 - T 2 )(3 2 At 2 d t . 2 d 2 Xl U 2 . (28) 

If n 2 = </>, then it follows from Eq. (JTHJ) that 

3T 3 /3 2 At 2 d t2 d 2 Xl U 2 + T 2 Atd 2 2 <j> = -r 2 (3r 3 - r|)/3 2 2 At 3 ^^. (29) 

Combining above expression and Eq. (|25[) . we modify II40 as 

n 40 = J [{d^if + 6/3 2 a n 2 (a n 1 ) 2 + 4/3 3 a n 3 ^ni + a 40 ] #, (so) 

where A 40 = I^I^lUK^ an d Eq. (27J> becomes 



9 t4 + a 4 ^ i n 4 = 0, (31) 

with «4 = Ai 3 T4/34. 

Therefore, when n 2 = 0, combining Eqs. (I15[) . (fT5|) . (f2"2"| and (|3ip . the fourth-order NPDE is exactly recovered to 
order 0(e 4 ) 

dt</> + dMt) + a 2 9 2 n 2 (0) + a 3 a 3 n 3 (0) + a 4 d 4 n 4 (0) = o, (32) 

with a 2 = Atr 2 f3 2 , a 3 = At 2 r 3 f3 3 , a 4 = At 3 r 4 /3 4 . 

Remark 2. If a 2 = or II 2 = 0, then we take j3 2 = which leads to Aio=0. For this case, the modification of H40 
is not needed, and the fourth-order NPDE is exactly recovered. 

Remark 3. The present model with a D1Q5 lattice can be used to simulate the fourth-order NPDEs with the form 
as Eq. ([T]) which contains some Kuramoto-Sivashinsky-type equations. Recently, following the idea in Ref. [20j two 
similar LBGK models with order 0(e 4 ) were given in Refs. [4!| and one for solving a class of three-order NPDEs 
which were first solved by the LBGK model in Ref. [2(|, and the other for the generalized Kuramoto-Sivashinsky 
equation. It can be easily found that the models in Refs. [4!| and [5(| do not satisfy the moments conditions ^ for 
m = 4. In fact, the so-called amending- function [ii|[5(| with order 0(At 2 ) implies that there is a stronger assumption 
on the nonlinear terms in EDF of the models. In addition, no comparisons were given to show the higher-order LBGK 
model in Ref. [49( superior to some lower-order LBGK models. 

E. Recovery of the Fifth-Order NPDE 
Summing Eq. $T2§ over j, and using Eqs. ©, ©, (fI5]l. (fT7| . and ([HI), we obtain 

d H <p + ir 3 At 2 d t3 d 2 Xl {f3 2 W 2 ) + 2T 2 Atd t2 d t3 cb + 4T 4 At 3 d t2 d 3 Xl (fJ 3 n 3 ) + r 5 At 4 ^ = . (33) 

j 

Using the similar procedure as Eq. (f25|) . we can obtain 

<^ U -- ( 34 ) 

j 
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with 



n 50 = J [(9 ni) 5 + io/3 2 ^n 2 (5 n!) 3 + 10 p^u^Ux) 2 + sp^n^iix] (35) 

and Eq. ([55)) becomes 

9 t5 + 3r 3 /3 2 Ai 2 flfefl^IIa + 2T 2 Atd t2 d t3 ct> + A Ti p 3 At 3 d t2 d 3 x n 3 + a 5 d 5 x U 5 = 0, (36) 



with «5 = Ai 4 T5/35. 

If II2 = II3 = </>, we can modify II50 as 



n 50 = / [(^no 5 + io/%fyn 2 (fyiii) 3 + lo^i^no 2 + 5/3 4 a n 4 a n 1 + a 50 ] d<p, (37) 



Where A 50 = ( 4 ^r 4 +3r|-2r|r 3 )/3 2 fe ^ and Eq ^ becomes 



9^0 + ^^5=0. (38) 
Combining Eqs. (fT5|). (fI5]l. d3TJ) and ([55]). the fifth-order NPDE is exactly recovered to order 0(e 5 ) 

5 

&0 + 3JIi(0) + ^a fc 3*n fc (0) =0, (39) 

with a fc = At k - 1 T k l3 k , k = 2,...,5. 

Remark 4- If a 2 = or 013 = 0, then j3 2 = or (3 3 = which leads to A50 = 0. For this case, the modification of 
II50 is not needed, and the fifth-order NPDE is exactly recovered. The present model with a D1Q6 or D1Q7 lattice 
can be used to simulate the fifth-order NPDE (JXJ) (m = 5) which contains some Kawahara-like equations. 

F. Recovery of the Sixth-Order NPDE 
Summing Eq. ([13]) over j, and using Eqs. ([3]), (0, ([151), El); (®, and ([25]), we obtain 

dt^ + 3T 3 At 2 d t4 d 2 Xi ((3 2 U 2 ) + 2T 2 Atd t2 d u ct> + 6r 4 At 3 d 2 2 d 2 Xi (f3 2 IL 2 ) + r 3 At 2 d 3 </> + r 2 Atd 2 J 

+4r 4 At 3 d t3 d 3 Xl (/W + hT b At A d t2 d A Xl (/3 4 n 4 ) + r 6 At 5 D %ff = °- (4°) 

j 

Using the similar procedure as Eq. (|25p . we can obtain 

£z^./j 0) =/36C n 6> (41) 

j 

with 

n 60 = / [(^n^ 6 + i5f3 2 dMdM 4 + 2o/3 3 ^n 3 (a n 1 ) 3 + 15 p^u^Ux) 2 + Q^d^d^x] d<p, (42) 

and Eq. (fH)]) becomes 

a 4e + 3T 3 /? 2 At 2 d t4 9^112 + 2T 2 Atd t2 d ti( j> + 6r^ 2 At 3 d 2 2 dlU 2 + r 3 At 2 d 3 J + r 2 Atd 2 J 

+4T i p 3 At 3 d t3 d 3 Xl Tl 3 + 5r 5/ 3 4 Ai 4 d 42 d^IFt + a 6 C n 6 = 0, (43) 



with a@ = At 5 TQ/3e. 

If II 2 = II 3 = II 4 = <fi, we can modify Ilgo as 



n 60 = f [(^n x ) 6 + ls/^n^nx) 4 + 2o/? 3 a n 3 (a n 1 ) 3 + ls&^i^nx) 2 + sfcd^d^x + a 60 ] (44) 
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Where A 60 = + + ^ ^ Eq becQmes 



d t6 <t> + a 6 d*Jle = 0. (45) 
Combining Eqs. (JUJ), JTHJ), d22j) , (JHTJ) , ([38]), and (45|) the sixth-order NPDE is exactly recovered to order 0(e 6 ) 



d t c/> + dxllx {<£) + Oikd k Jlk (0) = 0, 



(46) 



k=2 



with a k = At^Tkfa, k = 2, . . . , 6. 

Remark 5. If «2 = and H3 = </>, then the sixth-order NPDE is exactly recovered. The present model with a D1Q7 
lattice can be used to simulate six-order NPDE |T]) (m = 6) which also contains some Kawahara-like equations. 

III. EQUILIBRIUM DISTRIBUTION FUNCTIONS AND AUXILIARY MOMENTS 

For a given NPDE of order m, Hko can be obtained from the related formula in above section (2 < k < to), and 
the number of discrete velocity directions is at least equal to m + 1. So we can use a D1Q5 LBGK model to solve the 
NPDE of order less than or equal to 4, and a D1Q7 one to solve the NPDE of order less than or equal to 6. A D1Q4 
or D1Q6 one without the rest velocity can also be used to the NPDE of order 3 or 5. In this section we give only the 
EDFs and AMs for the LBGK models with D1Q5 and D1Q7 lattice, respectively. 

A. Equilibrium Distribution Functions for LBGK Model with D1Q5 Lattice 

Denoting n = <f>, Si = 111 /c, Sfc = (Hko + Pk^-k)/ c k , k = 2, 3, 4, the moments conditions ([3]) for m — 4 are rewriten 



E e "/r = n fe ,fc = o,...,4 

where {eo, ei, e2, e^, ca} = {0, 1, — 1, 2, —2}. Let 

ri = [no,n 1 ,...,n 4 ] T ,f cc '^[/ cq ,/ 1 cq , 

From Eq. (|4"7]) , we have 



eqiT 



,/4 q ] 



(47) 



(48) 



where 



It is easy to find the inverse of M5 



M, 



Mr 1 



1 

24 



M 5 f cq = 


n. 






"11 1 


1 


1 " 




1 -1 


2 


-2 




1 1 


4 


4 


5 


1 -1 


8 


-8 




.0 1 1 


16 


16. 




" 24 


-30 





6 


16 


16 


-4 


-4 


-16 


16 


4 


-4 


-2 


-1 


2 


1 


2 


-1 


-2 


1 



(49) 
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thus 

f cq = M^IL (50) 
Therefore, the EDFs of the LBGK model with D1Q5 lattice can be obtained as follows 

/ cq = [4^ - 5n 2 + n 4 ] /4, 

A cq = [4(fii + n 2 ) - n 3 - n 4 ] /6, 

/ 2 cq = [4(-n 1 + n 2 ) + n 3 -n 4 ] /e, 

/ 3 cq = [-2(n 1 -n 3 )-n 2 + n 4 ] /24, 

/ 4 cq = [2(n 1 -n 3 )-n 2 + n 4 ] /24,. (51) 



B. Equilibrium Distribution Functions for LBGK with D1Q7 Lattice 

6, the moments conditions ^ for m = 6 

(52) 



Similarly, denoting n = <p, U x = Ux/c, U k = (n k0 + /3 fc n fc )/c fe , k = 2 
are rewriten as 



^e 3 fc /; q = fl fc ,fc = 0,...,6 



where {e , ei, e 2 , e 3 , e 4 , e 5 , e 6 } = {0, 1,-1, 2, -2, 3, -3}. Let 

n = [n ,n 1 ,...,n 6 r,f eq = [/ eq ,/ 1 eq ,. 

From Eq. (|52|) . we have 

M 7 f eq = n, 

where 



M 7 = 



fCqiT 
> /6 J • 



It is easy to find the inverse of M7 



Mf 1 = — 
7 720 





" 1 


1 1 


1 


1 1 


1 " 











1 -1 


2 


-2 3 


-3 











1 1 


4 


4 9 


9 











1 -1 


8 


-8 27 


-27 











1 1 


1G 


16 81 


81 











1 -1 


32 - 


32 243 


-243 











1 1 


G4 


64 729 


729 _ 






720 





-980 





280 





-20" 





540 


540 


-195 - 


195 


15 


15 





-540 


540 


195 - 


195 - 


15 


15 





-108 


-54 


120 


60 - 


12 


-6 





108 


-54 


-120 


60 


12 


-6 





12 


4 


-15 


-5 


3 


1 





-12 


4 


15 


-5 - 


-3 


1 _ 






= M 7 


1 rl. 









thus 



Therefore, the EDFs of the LBGK model with D1Q7 lattice can be obtained as follows 

/ cq = [360 - 49fl 2 + 14S4 - fl 6 ] /36, 

A oq = [36(ni + n 2 )-i3(n3 + n 4 ) + n 5 + n 6 ]/48, 
[36(-fli + n 2 ) + i3(n 3 - n 4 ) - n 5 + n 6 ] /48, 
[-1811! - 9n 2 + 2on 3 + ion 4 - 2n 5 - n 6 ] /120, 
[i8ni - 9n 2 - 2on 3 + ion 4 + 2n 5 - n 6 ] /120, 



J2 
J3 



J4 

/ 5 oq = [1211! + 4n 2 - i5n 3 - 5n 4 + 3n 5 + n 6 ] /720, 
/ 6 oq = [-i2ni + 4n 2 + i5n 3 -5n 4 -3n 5 + n 6 ] /720, 



(53) 
(54) 



(55) 



(56) 
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C. Auxiliary Moments 

For a given NPDE of order m, 11^(1 < k < m) are known, and Ilfeo(2 < k < m) can be obtained by using the 
formula in Sec. II., then we can obtain the EDFs of related LBGK model with suitable lattice, and the LBGK model 
is derived. For example, the EDFs of D1Q5 or D1Q7 LBGK model can be computed by Eq. (j5"Tj) , or Eq. ([55|) . In 
simulations in the present work, two classes of concrete NPDEs are considered, and the related AMs of these NPDEs 
are given below. 

(1) Sixth-order NPDEs 



u t + auu x + (3u n u x + 2J ctkd^u = 0, 

fc=2 



(57) 



where a, (3, n and a k {2 < k < 6) are constants, and n > 0. Equation (|57|) contains many well-known equations, such 
as the KdV-Burgers equation, (m)KdV equation, Kuramoto-Sivashinsky equation, Kawahara equation, and some of 
their extensions. 

We can use a LBGK model with D1Q7 lattice to solve Eq. ([57]) . and the first seven moments are needed. Since 
III = §u 2 + -^-^u n+1 . = u, k > 2, we have =a« + Pu n , Tl k — 1, k > 2. Equation (|5T[) can be exactly recovered 
by using the LBGK model proposed above, and the exact expressions of AM functions nfeo(2 < k < 6) can be obtained 
as follows by using the formula in Sec. II. 



where 



n, 



tin 



n 20 = J (K) du, 
n 30 = J [(Kf + 
n 40 = J [(n;) 4 + 
n 50 = / [(n;) 5 ^ 
[(n;) 5 + 
[(n;) 6 + 
[(ni) 6 + 



30 2 K n 'i] du = J (ni) 3 dw + 3/? 2 ni, 

6^n^(ni) 2 + 4#,n£ni + A 4a ] du = [(ntfdu + 6/? 2 n 20 + 4/3^ + a 4(]U , 
io/? 2 n 2 (ni) 3 + lo&n^n;) 2 + s&n^ni + a 50 ] du 

lO^Cn'i) 3 ] du + 10/3 3 n 20 + 5/3 4 IIi + A 50 u, 

i5/? 2 n' 2 (ni) 4 + 2o/3 3 n^(n' 1 ) 3 + ls&n^n;) 2 + e^n^ni + a 60 ] du 

15/3 2 (ni) 4 + 20/3 3 (n' 1 ) 3 ] + 15/3 4 n 20 + 6/35IL + A 60 w, 



a fc = At K - L T k f3 k ,2< k < 6, 

_ t 2 (3t 3 -t 2 )/3 2 2 
^i 4 n — , 



T4 



A 



50 



(4r 2 r 4 + 3r| - 2r 2 r 3 )/3 2 /3 3 



7-5 



(5r 2 r 5 + 3r 3 r 4 - 2r|r 4 )/3 2 /3 4 + (4t 4 - t 2 t 3 )t 3 /3 2 + (t 2 t 3 - 6r 4 )r 2 /3. 3 

T6 



-460 = 



and it is easily obtained that 

k 



f{U[) k du = J2 I C 3 k (au) k -i (f3u n y du = u 



E 

3=0 



jn + k + 1 — j 



C£(au)*-'(/3u n ) 



,2 < fc < 6. 



(58) 



(59) 



(60) 
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For k — 2, 3, 4 we have 

(iLj 2 ^ = 

J (U[fdu = 

(ni) 4 du = 

(2) Third-order NPDEs 



i(cm) 2 



n + 2 



(au)(/3u n ) + 



1 



2n + 1 



08« 



n\2 



l(a U ) 3 + -JL (au ) 2 (/3 u ") + + -^—{Pu n f 

4 n + 3 2n+2 in + 1 



C9t» 



n\4 



(61) 



u t + aMj + f3u n u x + a 2 (u p ) xx + a 3 (u q ) xxx = 0, 



(62) 



where a,/?, a 2 , 03, n, p and <? are constants, and n > 0, p > 1, and q > 1. Eq. (|62p contains also many well-known 
equations, such as the KdV-Burgers equation, (m)KdV equation, and K(p, g) equation. 

We can use a LBGK model with D1Q4 or D1Q5 lattice to solve Eq. (|62|) . and if the C-E expansion to order 0(e 3 ) is 
used, we need only to give the first four moments, while for the C-E expansion to order 0(e 4 ), the first five moments 
are needed. Since 111 — ^u 2 + -J^u n+1 , II 2 = u p and n 3 = u q , we can exactly recover Eq. (|6"2"| by using the LBGK 
model proposed above, and the exact expressions of AM functions n 2 o, n 30 and LT40 can be obtained as follows by 
using the formula in Sec. II. 



n 20 
n 30 



n 



40 



(KYdu, 



[(ILJ 3 +3/3 2 n^ni] du = J [(U^) 3 + 3(3 2 pu p - 1 {au + f3u n )] du 
u (3u r > 



(U\) 3 du + 3(3 2 pu p 



p + 1 p+n 



J [(ilj 4 + 6/3 2 n 2 (n' x f + W 3 n[ + A iQ ] du 



p + n + 1 p + 2n 



n 1/ 



where A 40 , /(ILJ 2 ^, f(U[) 3 du, and /(ni) 4 dw can be computed by Eqs. ((59]) and lj6l]). 



IV. SIMULATION RESULTS 



q+1 q+n 



(63) 



To test the LBGK model proposed above, three classes of NPDEs with exact solutions are simulated, including four 
Kuramoto-Sivashinsky-like equations, three Kawahara-like equations and two KdV-like equations. In all simulations, if 
not specified, we use the nonequilibrium extrapolation scheme proposed by Guo et al. [5l| to treat the exact boundary 
condition, and the initial and boundary conditions of the test problems with analytical solutions are determined by 
their analytical solutions. The D1Q5 and D1Q7 LBGK models are used to simulate the test problems. The following 
global relative error is used to measure the accuracy: 



E = 



Ej\^,t)-4>*^j,t)\ 
E,- *)l 



(64) 



where cj> and <jf are the numerical solution and analytical one, respectively, and the summation is taken over all grid 
points. 

The first four test problems are the fourth-order Kuramoto-Sivashinsky-type equations, and three of them were 
simulated by the LBGK model in Ref. [50| . We use the D1Q5 LBGK model to simulate them and compare the 
proposed model with that in Ref. [5(| • 

Example 4-1- The Kuramoto-Sivashinsky equation [42| 



Uf -\- 1L1L X -\- U xx ~|- U xxxx 0, 



(65) 
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TABLE I: Comparison of global relative errors at different times [ (1) t = 1; (2) t = 2; (3) t = 3; (4) t = 4 ]. 



(c, T opt ) 



(10,5.99) 



Present Model 
(100,1. 



(1000, 1.27) 



(10,3.346) 



Model in Ref. [50] 
(100, 1.998) 



(1000,1.2705) 



(1) 
(2) 
(3) 
(4) 



9.6476 x 10 
1.2962 x 10 
1.7247 x 10 
2.2122 x 10 



2.9926 x 10 
4.2992 x 10 
5.5078 x 10 
6.8650 x 10 



6.0570 x 10 _ " 
8.7427 x 10~ 4 
1.1185 x 10~ 3 
1.3738 x 10~ 3 



2.6571 x 10 
3.7065 x 10 
5.0615 x 10 
8.0887 x 10 



3.2655 x 10 
5.3215 x 10 
7.1611 x 10 
8.9284 x 10 



6.5581 x 10 
1.1121 x 10 
1.5426 x 10 
1.9441 x 10 




10 



10 20 



FIG. 1: (Color online) The shock profile wave propagation of the KS equation (1651 1, 
boundary condition in [-30, 30], Ax = 0.1, At = 0.0001. 



5, k 



19 ' X ° 



-12. Exact 



with the exact solution 



u(x, i) 



-9tanh(fc(x - bt - x )) + 11 tanh 3 (fc(a; - bt - x ))) 



(66) 



where 6, fc, xq are parameters. 

In simulations, we set b = 5,k = \ \J^, and xq = —12 as in Ref. [5fj| for comparison. The simulation is conducted 
in [-30,30] with Ax = 0.1, At = 0.01,0.001, and 0.0001, corresponding to c = 10, 100, and 1000, respectively. The 
errors are listed in Table I for different times, where r opt is the optimal one corresponding to the minimal error. We 
also present the regular shock profile wave propagation for Eq. ([65]) with Eq. (|66]) in Fig. 1 . From the table it can be 
seen that the errors of our model are smaller than those of the model in Ref. [50( , and the accuracy of our model for 
At = 0.01 is much better than that of the model in Ref. [soj . When At is small enough, the effect of truncated errors 
of the model in Ref. can be ignored, thus the difference between the present model and that in Ref. (5(| is less. 

Example 4-2. The Kuramoto-Sivashinsky equation [42| 



with the exact solution 



u(x, t) = b + 



15 



19V19 



(—3 tanh(&(x — bt — xq)) + tanh 3 (A;(a: — bt — xq))) , 



(67) 



(68) 
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TABLE II: Comparison of global relative errors at different times [(!)*= 6; (2) t = 8; (3) t = 10; (4) t = 12 ]. 



Present Model Model in Ref. [50J 
(100,2.076) (1000,1.277) (10,4.0) (100,2.063) (1000, 1.277) 



(c,T opt ) (10,4.569) 



(1) 
(2) 
(3) 
(4) 



2.8486 x 10~ J 
3.1775 x 10" 5 
3.3937 x 10~ 5 
3.4934 x 10" 5 



1.5343 x 10 
1.6534 x 10 
1.7189 x 10 
1.8264 x 10 



2.7448 x 10" 
2.9535 x 10" 
3.0741 x 10" 
4.2625 x 10" 



1.2313 x 10 
1.5818 x 10 
1.9018 x 10 
2.1886 x 10 



5.6085 x 10~ J 
6.9643 x 10" 5 
8.1373 x 10~ 5 
9.1494 x 10~ 5 



7.6750 x 10" 
9.3058 x 10" 
1.0640 x 10" 
1.1661 x 10" 




FIG. 2: (Color online) The shock profile wave propagation of the KS equation (|67[) . b 
condition in [-80, 80], Ax = 0.1, At = 0.01. 



5,fc 



i 

2-/19 



,X 



-25. Exact boundary 



where b, k, xq are parameters. 

In simulations, we set b — 5, k — 9 Jr^ , and xq — —25 as in Ref. [501 ] for comparison. The simulation is conducted 
in [-50, 50] with Ax = 0.1, At = 0.01, 0.001 and 0.0001. The errors are listed in Table II. for different times, and the 



regular shock profile wave propagation for Eq. (|67[) is shown in Fig . 2. 
of our model are much smaller than those of the model in Ref. [50|. 
Example 4-3. The generalized Kuramoto-Sivashinsky equation [421 ] 



From the table it can be seen that the errors 



0, 



with the exact solution 



u(x, t) — b + 9 — 15 (tanh(fc(a; — bt — xq)) + tanh 2 (fc(a; — bt — Xq)) — tanh 3 (fe(a; — bt — Xq))) 



(69) 



(70) 



where a, b, k, xq are parameters. 

In simulations, we set a = 4, 6 = 6, k — |, and Xq = — 10 as in Ref. [5(| for comparison. The simulation is performed 
in [—30,30] with Ax = 0.1, At — 0.01,0.001 and 0.0001. Table III gives the errors of numerical solution at different 
times. We also present the solitary wave propagation for Eq. (f6"9")l is shown in Fig. 3. From the table it can be seen 
that the errors of our model for the smallest At are larger than those of the model in Ref. I50l , while the accuracy 
and stability of the present model for larger At are better than those of the model in Ref. [5fJ] . It is noted that the 
accuracy of both models for Eq. (I69|) is much larger than for Eqs . ([65| ) and (|67| . 

Example 4-4- The generalized Kuramoto-Sivashinsky equation [5 21 ] 



u t + 3u u x + au z 



bu r , 



0, 



(71) 
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TABLE III: Comparison of global relative errors at different times [ (1) t = 1; (2) t = 2; (3) t = 3; (4) t = 4; 
scheme is divergent ]. 



means that the 





(c,T opt ) (10,7.082) 




Present Model 
(100, 9.89) 


(1000, 1.267) 


(10,3.47) 


Model in Ref. [50] 
(100, 1.975) 


(1000, 1.267486) 


(1) 


4.1701 x 10" 


i 


5.2017 x 10~ 2 


5.1020 x 10 -2 


9.7859 x 10 _i 


1.3802 x 10 _i 


2.6054 x 10~ 2 


(2) 


1.2376 x 10" 





6.9440 x 10~ 2 


5.6700 x 10~ 2 




1.4077 x 10 _1 


2.8329 x 10~ 2 


(3) 


2.5757 x 10" 





9.7967 x 10~ 2 


5.1337 x 10~ 2 




1.7050 x 10" 1 


2.6802 x 10~ 2 


(4) 


3.4682 x 10" 





1.6776 x 10" 1 


6.5639 x 10" 2 




3.1488 x 10 _1 


3.5225 x 10" 2 





FIG. 3: (Color online) The solitary wave propagation of the GKS equation 
condition in [-30, 30], Ax = 0.1, At = 0.0001. 



b = 6, a = 4, k = §, a;o = — 10. Exact boundary 



with the exact solution 



u(x, t) 



V3b 
V2 



tanh 



VSb ( 29b 3 \ C 



b 

6' 



(72) 



where a,b,C are constants. 

In simulations, we set a = 1, b = 1, and C = 1. The simulation is performed in [—30,30] with xq = 0, Ax = 
0.1, At = 0.01,0.001 and 0.0001, and both D1Q5 and D1Q7 LBGK models are used. Table IV gives the errors of 
numerical solution at different times, and the regular shock profile wave propagation for Eq. (|71[) is shown in Fig. 4. 
From the table It found that the numerical solutions are agree well with the analytic ones, and the accuracy of Dlq7 
model is much better than that of D1Q5 one for smaller At. 

The next three test problems are the fifth-order Kawahara-like equations. We use the D1Q7 LBGK model to 
simulate them. Since the first six constraints on moments are enough for exactly recovering the fifth-order NPDEs in 
C-E expansion, it is interesting to compare the present exact model with C-E expansion to order 0(e e ) with one to 
order 0(e 5 ). We denote scheme 1 and scheme 2 for the model of order 0(e 6 ) and that of order 0(e 5 ), respectively. 
For simplification, we only take Tig = in Eq. (|56[) for scheme 2 in simulations. 

Example 4-5. The Kawahara equation [53| 



u t + auu x + (3u xxx + ju xxxxx = 0, 



(73) 
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TABLE IV: Comparison of global relative errors at different times [ (1) t = 1; (2) t = 2; (3) t = 3; (4) t = 4 ]. 





(c,T op t) (10,3.32) 




D1Q5 Model 
(100,2.0) 


(1000, 1.27) 


(10,4.14) 




D1Q7 Model 
(100,2.31) 


(1000, 1.40) 


(1) 


1.4921 x 10" 


-a 


3.7819 x 10" 4 


7.5629 x 10~ s 


1.3234 x 10" 


3 


8.7735 x 10~ a 


4.2032 x 10"° 


(2) 


3.1612 x 10" 


-3 


7.8311 x 10~ 4 


1.5509 x 10" 4 


2.6053 x 10" 


-3 


1.5272 x 10~ 4 


6.9690 x 10~ 6 


(3) 


5.0988 x 10" 


-3 


1.2215 x 10 -3 


2.4006 x 10~ 4 


4.1570 x 10" 


-3 


2.2473 x 10~ 4 


1.0021 x 10~ 5 


(4) 


7.2939 x 10" 


-3 


1.6930 x 10" 3 


3.3080 x 10~ 4 


6.0013 x 10" 


-3 


3.0868 x 10~ 4 


1.3604 x 10~ 5 




with the exact solution 

u(x,t) = Asech 4 ^^ - Ct)], (74) 

where a, f3, 7 are constants, and A = -jgH^, B = l^j-JL., C = -fff^. 

In simulations, we set a = (3 = —7 = 1. The simulation is conducted in [—30,30] with Ax = 0.1, At = 0.01,0.001 
and 0.0001. Table V. gives the errors of numerical solution for different times. From the table it is found that there 
is little difference between accuracy of two schemes, which implies that the accuracy of higher-order model may not 
be better than that of lower-order one. 



TABLE V: Comparison of global relative errors at different times [(!)*= 1, (2) t = 2; (3) t = 3; (4) t = 4 ]. 





(C> Topt) 


(10,3.37) 




Scheme 1 
(100,2.53) 




(1000,2.02) 


(10,3.35) 




Scheme 2 
(100,2.55) 




(1000, 2.04) 


(1) 




6.0101 x 10" 


3 


2.6928 x 10" 


y 


1.5361 x 10~ 3 


5.9364 x 10- 


3 


2.7372 x 10" 


3 


1.5750 x 10" 3 


(2) 




1.0877 x 10" 


2 


5.2590 x 10" 


3 


3.0032 x 10 -3 


1.0698 x 10" 


2 


5.3477 x 10" 


3 


3.0827 x 10~ 3 


(3) 




1.5605 x 10" 


-2 


7.5403 x 10" 


3 


4.2350 x 10~ 3 


1.5369 x 10" 


2 


7.6869 x 10" 


3 


4.3599 x 10" 3 


(4) 




2.0197 x 10" 


-2 


9.4960 x 10" 


3 


5.3288 x 10~ 3 


2.0035 x 10" 


2 


9.7447 x 10" 


3 


5.4829 x 10 -3 



15 





TABLE VI: Comparison 


of global relative 


errors at different times 


[ (1) t = 1, (2)t 


= 2; (3) t = 3; 


(4)i 


= 4]. 








Scheme 1 






Scheme 2 








(c, Topt) (10, 4.54) 




(100,2.53) 


(1000,2.04) 


(10,4.04) 


(100,2.56) 




(1000, 2.04) 


(1) 


1.9295 x 1(T 


-2 


7.1698 x 10" 3 


4.4255 x 10 -3 


1.6850 x 10'* 


7.3109 x 10" 


a 


4.4254 x 10~ J 


(2) 


3.8260 x 10" 


-2 


1.3214 x 10"" 2 


7.7342 x 10~ 3 


3.2649 x 10~ 2 


1.3549 x 10" 


-2 


7.7745 x 10~ 3 


(3) 


5.7488 x 10" 


-2 


1.7575 x 10" 2 


1.0046 x 10" 2 


4.7022 x 10~ 2 


1.8137 x 10" 


-2 


1.0211 x 10~ 2 


(4) 


7.4409 x 10" 


-2 


2.1032 x 10~ 2 


1.1886 x 10" 2 


6.0945 x 10~ 2 


2.2855 x 10" 


2 


1.2424 x 10~ 2 



TABLE VII: Comparison of global relative errors at different times [ (1) t = 1, (2) t = 2; (3) t = 3; (4) t = 4 ]. 





(C, T opt ) 


(10,5.01) 




Scheme 1 
(100,3.31) 




(1000,2.28) 


(10,4.64) 




Scheme 2 
(100,2.95) 




(1000,2.18) 


(1) 




9.8169 x 10" 


-3 


4.6248 x 10" 


3 


2.0662 x 10~ 3 


9.0008 x 10" 


-3 


3.6976 x 10" 


-3 


1.8557 x 10~ 3 


(2) 




1.8335 x 10" 


-2 


8.9484 x 10" 


3 


4.0584 x 10~ 3 


1.6575 x 10" 


-2 


7.1666 x 10" 


■3 


3.6371 x 10~ 3 


(3) 




2.6841 x 10" 


-2 


1.3201 x 10" 


2 


6.1428 x 10~ 3 


2.3973 x 10" 


-2 


1.0468 x 10" 


-2 


5.2117 x 10~ 3 


(4) 




3.5872 x 10" 


-2 


1.8257 x 10" 


2 


7.9304 x 10" 3 


3.1522 x 10" 


-2 


1.4085 x 10" 


-2 


6.6671 x 10~ 3 



Example The modified Kawahara equation [54| 

Ut + flu u x + bu xxx ku xxxxx = 0, (75) 

with the exact solution 

u(x, t) = Asech 2 [B (x - Ct)] , (76) 
where a. b, k are constants, and A = — 7==, B = i\/ -X-.C = 

' ' ' VlOa*: ' 2 y 5k ' 25k 

In simulations, we set a = b = k = 1. The simulation is conducted in [—30, 30] with Aa; = 0.1, At = 0.01, 0.001 and 
0.0001. Table VI. gives the errors of numerical solution for different times. From the table it is also found that there 
is little difference between accuracy of two schemes. 

Example 4-7. The Korteweg-de Vries-Kawahara equation [iij 

~t~ Ullx ~\~ Ux "T" U XX x Uxxxxx — 0, (77) 

with the exact solution 

, . 105 , 4 / 1 . 205 A 

u(x,t) — sech — i=(x 1 — xn) ) , (78) 

169 \2Vl3 169 u 7 ' v ; 

where xo is a parameter. 

In simulations, we set xo — 20 as in Ref. Q. The simulation is conducted in [0, 200] with Aa; = 0.1, At — 0.01, 0.001 
and 0.0001. Table VII. gives the errors of numerical solution for different times. It can be found that the numerical 
solutions obtained by LBGK model are in good agreement with the analytic ones. From the table it is still found that 
there is little difference between accuracy of two schemes. 

The last two test problems are the third-order KdV-type equations. We use the D1Q5 LBGK model to simulate 
them. Similar to the above three test problems, we also compare the present exact model with C-E expansion to 
order 0(e 4 ) with one to order 0(e 3 ), and denote scheme 1 and scheme 2 for the model of order 0(e 4 ) and that of 
order 0(e 3 ), respectively. For simplification, we only take II4 = in Eq. (|51[) for scheme 2 in simulations. 

Example 4-8- The KdV Burgers equation [2(| 



utt + auu x - ju X x + Suxxx = 0, (79) 



with the exact solution 



where v = --^ and £ = fgy. 



u( M ) = 2£(l- [i + eMa _ gf)]2 j, (80) 
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TABLE VIII: Comparison of global relative errors for c = 1 at different times [ (1) Scheme 1; (2) Scheme 2; (3) Scheme in Ref. 
01- 



7~opt 




t = 10 t = 50 t = 150 




£ = 250 


t = 300 


0.97 
0.96 
0.968 


(1) 
(2) 
(3) 


2.4300 x 10~ B 4.2738 x 10" B 4.0518 x 10" 
4.9172 x 10~ 6 7.1775 x 10" 6 6.1615 x 10" 
1.0416 x 10" 5 1.8801 x 10" 5 1.7409 x 10" 


6 
■5 


3.4676 x 10" b 
5.2459 x 10" 6 
1.4877 x 10" 5 


3.3242 x 10"° 
4.9383 x 10~ 6 
1.3901 x 10" 5 






TABLE IX: Comparison of global relative errors for 


C 


= 10 at different times. 




Topt 




t = 1 t = 2 




t = 3 


t = 4 


37.77 
31.52 


Scheme 1 1.8629 x lO" 3 9.1678 x 10" 4 
Scheme 2 2.6472 x 10" 2 1.0947 x 10" 2 




9.2998 x 10" 4 
6.1720 x 10~ 3 


7.1608 x 10" 4 
3.3151 x 10 -3 



In simulations, we set a = 1, 7 = 9 x 10~ 4 , 6 = 2 x 10~ 5 , and the simulation is conducted in [—4, 4] with Ax — At — 
0.01 as in Refs. [2(| and [49[. Table VIII. gives the errors of numerical solution for different times. From the table it 
is found that the errors of schemes 1 and 2 are much smaller than those of the model in Ref. [4!| , and scheme 1 is 
better that scheme 2, but there is little difference between our schemes. The numerical results show that the model 
in Ref. [49| is not exact, even for the constraints on lower-order moments. 

Example 4-9. The K(n, n)-Burgers equation [55| 

u t + a(u n ) x + b{u n ) xxx + ku xx = 0, (81) 

with the exact solution 

u(x, t) = [A(l + tanh(Bx + Ct))]~~^ ,6>0,n>l,a<0 (82) 

where a, b, k, n are constants, and A = j^V—ab, B = C = "^"^ ■ 

In simulations, we set a = —1, b = 1, k = — 1, n = 2. The simulation is conducted in [— 1, 1] with Ax = 0.01, At = 
0.001. Table IX. gives the errors of numerical solution for different times. From the table it is found that the errors of 
schemes 1 are much smaller than those of scheme 2, which implies that for this problem the accuracy of higher-order 
model is much better than that of lower-order one. 



V. CONCLUSION 

In the present work, we have developed a unified LBGK model for ID higher-order NPDEs. Through C-E expansion 
a given NPDE can be exactly recovered to required order of small parameter e by choosing correct auxiliary moments. 
Unlike traditional numerical methods which solve for macroscopic variables, the model has the advantages of standard 
LBGK model, which are borrowed from kinetic theory, such as linearity of the convection operator in velocity space, 
simplicity and symmetry of scheme, ease in coding and intrinsical parallelism [3[. Detailed numerical tests of the 
proposed model are carried out for different types of NPDEs, including the Kuramoto-Sivashinsky type equations, 
Kawahara type equations, and KdV type equations. It is found that the simulation results agree well with the 
analytical and numerical solutions reported in previous studies, which shows that the LBM has potentials in simulating 
higher-order NPDEs. However, perhaps due to the effects of nonlinearity and higher order differentials, the LBGK 
model for solving higher order NPDEs is sensitive to the key parameters, such as Ax, At and r, and it does not seems 
so efficient as that for solving lower order ones, such as that for NCDEs 39]. 

Note that the proposed model can be directly applied to derive the LBGK model for high-order NPDEs in higher 
dimensional space by treating moments as tensors, and the LBGK model for NPDEs with order larger than six can 
be easily derived by using the idea in this paper. Nevertheless, some important issues, such as how to improve the 
accuracy and stability of the LB models need further studies. 
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